      SUBROUTINE STEP (J1,J2,DT,ALPH,ROB)
C--
C--   SUBROUTINE STEP (J1,J2,DT,ALPH,ROB)
C--
C--   Purpose: perform one time step starting from u(1) 
C--            and using the 4th order Runge-Kutta scheme:
C--    h(1)  = F(u(1))
C--    u(2)  = u(1) + (DT/2)*h(1)
C--    h(2)  = F(u(2))
C--    u(3)  = u(1) + (DT/2)*h(2)
C--    h(3)  = F(u(1))
C--    u(4)  = u(1) + DT*h(3)
C--    h(4)  = F(u(4))
C--    G     = (1/6)*(h(1)+2*h(2)+2*h(3)+h(4))
C--   unew  = u(1) + DT*G 
C--
C--   In this scheme,  all tendency terms, including advection, diffusion 
C--   and gravity-wave terms, are treated explicitly with RK4.
C--
C--   **** This routine is for validation purpose. Nonlinear terms are 
C--   **** switched off, and the result is to be compared with Lorenz 4-cycle
C--   ****                                                      2012/11/29 DH
C--
C--   Input: 
C--   J1, J2     : dummy
C--   DT = time step (if DT < or = 0, tendencies are computed but 
C--                   no time stepping is performed)
C--   ALPH, ROB  : dummy
C-- 
C--   Modified common blocks : DYNSP1, DYNSP2
C--
C--                                           2012/11/13 Daisuke Hotta
      include "atparam.h"
      include "atparam1.h"

      include "com_dyncon0.h"
      include "com_hdifcon.h"

      include "com_dynvar.h"

      COMPLEX VORDT(MX,NX,KX), DIVDT(MX,NX,KX), TDT(MX,NX,KX),
     *        PSDT(MX,NX), TRDT(MX,NX,KX,NTR)
      COMPLEX VORH1(MX,NX,KX), DIVH1(MX,NX,KX), TH1(MX,NX,KX), 
     *        PSH1(MX,NX), TRH1(MX,NX,KX,NTR)                  
      COMPLEX VORH2(MX,NX,KX), DIVH2(MX,NX,KX), TH2(MX,NX,KX), 
     *        PSH2(MX,NX), TRH2(MX,NX,KX,NTR)                  
      COMPLEX VORH3(MX,NX,KX), DIVH3(MX,NX,KX), TH3(MX,NX,KX), 
     *        PSH3(MX,NX), TRH3(MX,NX,KX,NTR)                  
      COMPLEX VORH4(MX,NX,KX), DIVH4(MX,NX,KX), TH4(MX,NX,KX), 
     *        PSH4(MX,NX), TRH4(MX,NX,KX,NTR)                  

      iitest=0
      if(iitest.eq.1) print*, ' inside step'
      print *, 'DBG: delt, alpha' , DT, ALPH                   ! DH

C--   1. Computation of h1 (= F_exp(u(1)) ) and u(2)
      if (iitest.eq.1) print*,' call grtend'
c     CALL GRTEND      (VORH1,DIVH1,TH1,PSH1,TRH1,1,1)
      CALL CLEARTEND   (VORH1,DIVH1,TH1,PSH1,TRH1)
      CALL SPTEND            (DIVH1,TH1,PSH1,1)
      CALL HORDIF_ALL(1,VORH1,DIVH1,TH1,     TRH1)
      CALL TIMINT_ALL(DT*0.5, 2, VORH1, DIVH1, TH1, PSH1, TRH1)

C--   2. Computation of h2 (= F_exp(u(2)) ) and u(3)
      if (iitest.eq.1) print*,' call grtend'
c     CALL GRTEND      (VORH2,DIVH2,TH2,PSH2,TRH2,2,2)
      CALL CLEARTEND   (VORH2,DIVH2,TH2,PSH2,TRH2)
      CALL SPTEND            (DIVH2,TH2,PSH2,2)
      CALL HORDIF_ALL(2,VORH2,DIVH2,TH2,     TRH2)
      CALL TIMINT_ALL(DT*0.5, 3, VORH2, DIVH2, TH2, PSH2, TRH2)

C--   3. Computation of h3 (= F_exp(u(3)) ) and u(4)
      if (iitest.eq.1) print*,' call grtend'
c     CALL GRTEND      (VORH3,DIVH3,TH3,PSH3,TRH3,3,3)
      CALL CLEARTEND   (VORH3,DIVH3,TH3,PSH3,TRH3)
      CALL SPTEND            (DIVH3,TH3,PSH3,3)
      CALL HORDIF_ALL(3,VORH3,DIVH3,TH3,     TRH3)
      CALL TIMINT_ALL(    DT, 4, VORH3, DIVH3, TH3, PSH3, TRH3)

C--   4. Computation of h4 (= F_exp(u(4)) )
c     CALL GRTEND      (VORH4,DIVH4,TH4,PSH4,TRH4,4,4)
      CALL CLEARTEND   (VORH4,DIVH4,TH4,PSH4,TRH4)
      CALL SPTEND            (DIVH4,TH4,PSH4,4)
      CALL HORDIF_ALL(4,VORH4,DIVH4,TH4,     TRH4)

C--   5. Computation of **DT
      DO K=1,KX
         DO M=1,MXNX
            VORDT(M,1,K) = (1./6.)*(     VORH1(M,1,K) + VORH4(M,1,K) 
     *                             + 2.*(VORH2(M,1,K) + VORH3(M,1,K)))
            DIVDT(M,1,K) = (1./6.)*(     DIVH1(M,1,K) + DIVH4(M,1,K) 
     *                             + 2.*(DIVH2(M,1,K) + DIVH3(M,1,K)))
            TDT(M,1,K)   = (1./6.)*(       TH1(M,1,K) +   TH4(M,1,K) 
     *                             + 2.*(  TH2(M,1,K) +   TH3(M,1,K)))
            DO ITR=1,NTR
            TRDT(M,1,K,ITR)=(1./6.)*(TRH1(M,1,K,ITR)+TRH4(M,1,K,ITR) 
     *                         +2.*(TRH2(M,1,K,ITR)+TRH3(M,1,K,ITR)))
            END DO
         END DO
      END DO

      DO M=1,MXNX
      PSDT(M,1)=(1./6.)*(PSH1(M,1)+PSH4(M,1)+2.*(PSH2(M,1)+PSH3(M,1)))
      END DO

C--   7. Time integration

      IF (DT.LE.0.) RETURN

      if (iitest.eq.1) print*,' time integration'
      CALL TIMINT_ALL(    DT, 1, VORDT, DIVDT, TDT, PSDT, TRDT)

C--
      RETURN
      END   

      SUBROUTINE HORDIF (NLEV,FIELD,FDT,DMP)
C--
C--   Aux. subr. HORDIF (NLEV,FIELD,FDT,DMP)
C--   Purpose : Add horizontal diffusion tendency of FIELD 
C--             to spectral tendency FDT at NLEV levels
C--             using damping coefficients DMP
C--
      include "atparam.h"

      COMPLEX FIELD(MXNX,NLEV), FDT(MXNX,NLEV)
      REAL    DMP(MXNX)

      DO K=1,NLEV
        DO M=1,MXNX
          FDT(M,K)=FDT(M,K)-DMP(M)*FIELD(M,K)
        ENDDO
      ENDDO

      RETURN
      END

      SUBROUTINE HORDIF_ALL (J1,VORH,DIVH,TH,TRH)
C--
C--   Aux. subr. HORDIF_ALL (J1,VORH,DIVH,TH,TRH)
C--   Purpose : Apply horizontal diffusion to all variables
C--   Input:
C--     J1: time index
C--   Input/Output:
C--     VORH,DIVH, TH, TRH : tendencies to be modified

      include "atparam.h"
      include "atparam1.h"
      include "com_dyncon0.h"
      include "com_hdifcon.h"
      include "com_dynvar.h"

      INTEGER J1
      COMPLEX VORH(MX,NX,KX),DIVH(MX,NX,KX)
      COMPLEX TH(MX,NX,KX),TRH(MX,NX,KX,NTR)

      COMPLEX CTMP(MX,NX,KX)

C     1 Diffusion of wind and temperature
 
      CALL HORDIF (KX,VOR(1,1,1,J1),VORH,DMP)
      CALL HORDIF (KX,DIV(1,1,1,J1),DIVH,DMPD)

      DO K=1,KX
        DO M=1,MXNX
          CTMP(M,1,K)=T(M,1,K,J1)+TCORH(M,1)*TCORV(K)
        ENDDO
      ENDDO

      CALL HORDIF (KX,CTMP,TH,DMP)

C    2 Stratospheric diffusion and zonal wind damping

      SDRAG=1./(TDRS*3600.)
      DO N=1,NX
        VORH(1,N,1)=VORH(1,N,1)-SDRAG*VOR(1,N,1,J1)
        DIVH(1,N,1)=DIVH(1,N,1)-SDRAG*DIV(1,N,1,J1)
      ENDDO

      CALL HORDIF (1,VOR(1,1,1,J1), VORH,DMPS)
      CALL HORDIF (1,DIV(1,1,1,J1), DIVH,DMPS)
      CALL HORDIF (1,CTMP,TH ,DMPS)

C     3 Diffusion of tracers

      DO K=1,KX
        DO M=1,MXNX
          CTMP(M,1,K)=TR(M,1,K,J1,1)+QCORH(M,1)*QCORV(K)
        ENDDO
      ENDDO

      CALL HORDIF (KX,CTMP,TRH,DMPD)

      IF (NTR.GT.1) THEN
        DO ITR=2,NTR
          CALL HORDIF (KX,TR(1,1,1,J1,ITR),TRH(1,1,1,ITR),
     &                 DMP)
        ENDDO
      ENDIF
      END

      SUBROUTINE TIMINT (DT,NLEV,FOLD,FNEW,FDT)
C--
C--   Aux. subr. TIMINT (DT,NLEV,FOLD,FNEW,FDT)
C--   Purpose : Perform time integration of FOLD to FNEW at NLEV levels
C--             using tendency FDT
C--                          2012/11/13 Daisuke Hotta
      include "atparam.h"

      COMPLEX FOLD(MXNX,NLEV), FNEW(MXNX,NLEV), FDT(MXNX,NLEV)
      REAL DT
      INTEGER NLEV

      IF (IX.EQ.IY*4) THEN
        DO K=1,NLEV
          CALL TRUNCT (FDT(1,K))
        ENDDO
      ENDIF

      DO K=1,NLEV
        DO M=1,MXNX
          FNEW(M,K) = FOLD(M,K)+DT*FDT(M,K)
        ENDDO
      ENDDO

      RETURN
      END
C--  
      SUBROUTINE TIMINT_ALL(DT,J,VORH,DIVH,TH,PSH,TRH)
C--
C--   Aux. subr. TIMINT_ALL(DT,J,VORH,DIVH,TH,PSH,TRH)
C--   Purpose : Perform time integration for all variables
C--             using tendencies PSH,VORH,DIVH,TH,TRH
C--                          2012/11/13 Daisuke Hotta
C--
      include "atparam.h"
      include "atparam1.h"
      include "com_dynvar.h"
      
      COMPLEX VORH(MX,NX,KX), DIVH(MX,NX,KX), TH(MX,NX,KX),
     *        PSH(MX,NX), TRH(MX,NX,KX,NTR)
      REAL    DT
      INTEGER J
      
      CALL TIMINT (DT, 1,   PS(1,1,1),   PS(1,1,J), PSH)
      CALL TIMINT (DT,KX,VOR(1,1,1,1),VOR(1,1,1,J),VORH)
      CALL TIMINT (DT,KX,DIV(1,1,1,1),DIV(1,1,1,J),DIVH)
      CALL TIMINT (DT,KX,  T(1,1,1,1),  T(1,1,1,J),  TH)
      DO ITR=1, NTR
         CALL TIMINT(DT, KX, TR(1,1,1,1,ITR), TR(1,1,1,J,ITR),
     *               TRH(1,1,1,ITR))
      END DO
      
      RETURN
      END

      SUBROUTINE CGRATE (VOR,DIV,VORDT,DIVDT)
C--
C--   SUBROUTINE CGRATE (VOR,DIV,VORDT,DIVDT)
C--
C--   Purpose: Check growth rate of eddy kin. energy 
C--   Input  : VOR    = vorticity
C--            DIV    = divergence
C--            VORDT  = time derivative of VOR
C--            DIVDT  = time derivative of DIV
C--
      include "atparam.h"
      include "atparam1.h"

      COMPLEX VOR(MX,NX,KX), VORDT(MX,NX,KX), 
     &        DIV(MX,NX,KX), DIVDT(MX,NX,KX), TEMP(MX,NX)

      GRMAX=0.2/(86400.*2.)

      CDAMP=0.

      DO K=2,KX

        GRATE=0.
        RNORM=0.

        CALL INVLAP (VOR(1,1,K),TEMP)

        DO N=1,NX
         DO M=2,MX
           GRATE=GRATE-REAL(VORDT(M,N,K)*CONJG(TEMP(M,N)))
           RNORM=RNORM-REAL(  VOR(M,N,K)*CONJG(TEMP(M,N)))
         ENDDO
        ENDDO

        IF (GRATE.GT.GRMAX*RNORM) CDAMP =
     &      MAX(CDAMP,0.8*GRATE/RNORM)
C    &      MAX(CDAMP,(GRATE*GRATE)/(GRMAX*RNORM*RNORM))

      ENDDO

      IF (CDAMP.GT.0.) THEN

        print *, ' rot. wind damping enabled'

        DO K=1,KX
          DO N=1,NX
           DO M=2,MX
             VORDT(M,N,K)=VORDT(M,N,K)-CDAMP*VOR(M,N,K)
           ENDDO
          ENDDO
        ENDDO

      ENDIF


      CDAMP=0.

      DO K=2,KX

        GRATE=0.
        RNORM=0.

        CALL INVLAP (DIV(1,1,K),TEMP)

        DO N=1,NX
         DO M=2,MX
           GRATE=GRATE-REAL(DIVDT(M,N,K)*CONJG(TEMP(M,N)))
           RNORM=RNORM-REAL(  DIV(M,N,K)*CONJG(TEMP(M,N)))
         ENDDO
        ENDDO

        IF (GRATE.GT.GRMAX*RNORM) CDAMP =
     &      MAX(CDAMP,0.8*GRATE/RNORM)
C    &      MAX(CDAMP,(GRATE*GRATE)/(GRMAX*RNORM*RNORM))

      ENDDO

      IF (CDAMP.GT.0.) THEN

        print *, ' div. wind damping enabled'

        DO K=1,KX
          DO N=1,NX
           DO M=2,MX
             DIVDT(M,N,K)=DIVDT(M,N,K)-CDAMP*DIV(M,N,K)
           ENDDO
          ENDDO
        ENDDO

      ENDIF

C--
      RETURN
      END
C-- 
      SUBROUTINE CLEARTEND(VORH,DIVH,TH,PSH,TRH)
C--
C--   Aux. subr.  CLEARTEND(VORH,DIVH,TH,PSH,TRH)
C--   Purpose : Initialize tendencies (zero-clear)1
C--                          2012/11/28 Daisuke Hotta
C--
      include "atparam.h"
      include "atparam1.h"

      COMPLEX VORH(MX,NX,KX), DIVH(MX,NX,KX), TH(MX,NX,KX),
     *        PSH(MX,NX), TRH(MX,NX,KX,NTR)
      COMPLEX ZERO=(0.,0.)

      DO K=1,KX
        DO N=1,NX
          DO M=1,MX
            VORH(M,N,K)=ZERO
            DIVH(M,N,K)=ZERO
            TH  (M,N,K)=ZERO
          ENDDO
        ENDDO
      ENDDO

      DO ITR=1,NTR
        DO K=1,KX
          DO N=1,NX
            DO M=1,MX
              TRH (M,N,K,ITR)=ZERO
            ENDDO
          ENDDO
        ENDDO
      ENDDO

      DO N=1,NX
        DO M=1,MX
          PSH (M,N)=ZERO
        ENDDO
      ENDDO

      RETURN
      END
